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Extinction appears ubiquitously in many fields, including chemical reactions, pop- 
ulation biology, evolution, and epidemiology Even though extinction as a random 
process is a rare event, its occurrence is observed in large finite populations. Ex- 
tinction occurs when fluctuations due to random transitions act as an effective 
force which drives one or more components or species to vanish. Although there are 
many random paths to an extinct state, there is an optimal path that maximizes 
the probability to extinction. In this article, we show that the optimal path is asso- 
ciated with the dynamical systems idea of having maximum sensitive dependence 
to initial conditions. Using the equivalence between the sensitive dependence and 
the path to extinction, we show that the dynamical systems picture of extinction 
evolves naturally toward the optimal path in several stochastic models of epidemics. 
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1. Introduction 



Determining the conditions for epidemic extinction is an important public health 
problem. Global eradication of an infectious disease has rarely been achieved, but it 
continues to be a public health goal for polio (World Health Organization 20091) and 



many other diseases, including childhood diseases. More commonly, disease extinc- 
tion, or fade out, is local and may be followed by a reintroduction of the disease from 
other regions (lAnderson & May 1991 Grassly et al.\ 2005). Extinction may also 



occur for individual strains of a multistrain disease flMinayev & Ferguson 2009 1 



such as influenza or dengue fever. Since extinction occurs in finite populations, it 



depends critically on local community size ( Bartlett 1957 1960 Keeling & Grenfell 



1997 Conlan & Grenfell 20071. Moreover, it is important to know how parameters 



affect the chance of extinction for predicting the dynamics of outbreaks and for de- 



veloping control strategies to promote epidemic extinction ( Dykman et al. 2008 ) 



The determination of extinction risk is also of interest in related fields, such as evo- 
lution and ecology. For example, in the neutral theory of ecology, bio-diversity arises 
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from the interplay between the introduction and extinction of species (Banavar & 



Maritan 2009 Azaele et al. 2006). 



In general, extinction occurs in discrete, finite populations undergoing stochas- 
tic effects due to random transitions or perturbations. The origins of stochasticity 
may be internal to the system or may arise from the external environment. Small 
population size, low contact frequency for frequency-dependent transmission, com- 
petition for resources, and evolutionary pressure (de Castro & Bolker [2005]) , as 



well as heterogeneity in populations and transmission ( Lloyd et al. 
be determining factors for extinction to occur. 



2007), may all 



Extinction risk is affected by the nature and strength of the noise ( Melbourne 



& Hastings 


et al. 




200G 



2008), as well as other factors, including outbreak amplitude (Alonso 



and seasonal phase occurrence (Stone et al. 2007). For large popula- 



tions, the intensity of internal population noise is generally small. However, a rare, 
large fluctuation can occur with non-zero probability and the system may be able 
to reach the extinct state. For respiratory diseases such as SARS, super-spreading 
may account for these large stochastic fluctuations ( Lloyd-Smith et al. 2005 ) . Since 
the extinct state is absorbing due to effective stochastic forces, eventual extinction 



is guaranteed when there is no source of reintroduction (Bartlett 1949 Allen & 



Burgin 2000 Gardiner 2004 1 . However, because fade outs are usually rare events 



in large populations, typical time scales for extinction may be extremely long. 

Stochastic population models of finite populations, which include extinction 
processes, are effectively described using the master equation formalism. Stochas- 
tic master equations are commonly used in statistical physics when dealing with 
chemical reaction processes (Kubo 1963) and predict probabilities of rare events 
occurring at certain times. For many problems involving extinction in large popula- 
tions, if the probability distribution of the population is approximately stationary, 
the probability of extinction is a function that decreases exponentially with increas- 
ing population size. The exponent in this function scales as a deterministic quantity 
called the action (Kubo et al. 1973). It can be shown that a trajectory that brings 
the system to extinction is very likely to lie along a most probable path, called 
the optimal extinction trajectory or optimal path. It is a remarkable property that 
a deterministic quantity such as the action can predict the probability of extinc- 



tion, which is inherently a stochastic process ( Schwartz et al. 2009 Dykman et al 



2008). 



Locating the optimal path is desirable because the quantity of interest, the 
extinction rate, depends on the probability to traverse this path, and the effect 
of a control strategy on extinction rate can be determined by its effect on the 
optimal path (Dykman et al. 2008). By employing an optimal path formalism, 



we convert the stochastic problem to a mechanistic dynamical systems problem. In 
contrast to approaches based on diffusive processes that are valid only in the limit of 
large system sizes (Nasell 1999; Lind holmj 20081, this dynamical systems approach 
can give accurate estimates for the extinction time even for small populations if 
the action is sufficiently large. Additionally, unlike other methods that are used to 
estimate lifetimes, this approach enables one both to estimate lifetimes and to draw 
conclusions about the path taken to extinction. This more detailed understanding 



of how extinction occurs may lead to new stochastic control strategies (Dykman 



et al. 2008). 



In this article, we show that locating the optimal extinction trajectory can be 
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achieved naturally by evolving a dynamical system that converges to the optimal 
path. The method is based on computing finite-time Lyapunov exponents (FTLE), 



which have previously been used to find coherent structures in fluid flows (Haller 



2000 



2010 



2001] [20021 |Shadden et at] [20051 |Lekien et al\ [20071 |Branicki fc Wig gins 



The FTLE provides a measure of how sensitively the system's future be- 
havior depends on its current state. We argue that the system displays maximum 
sensitivity near the optimal extinction trajectory, which enables us to dynamically 
evolve toward the optimal escape trajectory using FTLE calculations. For several 
models of epidemics that contain internal or external noise sources, we illustrate 
the power of our method to locate optimal extinction trajectories. Although our 
examples are taken from infectious disease models, the approach is general and is 
applicable to any extinction process or escape process. 



2. Stochastic modeling in the large population limit 

To introduce our idea of dynamically constructing an optimal path to extinction in 
stochastic systems, we show its application to a stochastic Susceptible-Infectious- 
Recovered (SIR) epidemiological model. Details of the SIR model can be found in 
Sec. 1 of the electronic online supplement. Figure [l] shows the probability density of 
extinction prehistory in the SI plane. The probability density was numerically com- 
puted using 20, 000 stochastic SIR trajectories that ended in extinction. Trajectories 
are aligned at their extinction point. From the extinction point, the prehistory of 
each trajectory up to the last outbreak of infection is considered. Small fluctuations 
in the infectious population are not considered in identifying the last outbreak. In 
this way, we restrict the analysis to the interval between the last large outbreak 
of infection and the extinction point. The resulting (S,I) pairs of susceptible and 



infectious individuals are then binned and plotted in the SI plane (Eilers & Goe 



man 



2004 ) . The resulting discrete density has been color coded so that the brighter 
regions correspond to higher density of trajectories. The figure shows that, among 
all the paths that the stochastic system can take to reach the extinct state, there is 
one path that has the highest probability of occurring. This is the optimal path to 
extinction. One can see that the optimal path to extinction lies on the peak of the 
probability density of the extinction prehistory. It should be noted that extinction 



for the stochastic SIR model has been studied previously (Kamenev & Meerson 



2008) 



The optimal path can be obtained using methods of statistical physics. In fig- 
ure [TJ the numerical prediction of the entire optimal trajectory for the stochastic 
SIR system has been overlaid on the probability density of extinction prehistory 
that was found using stochastic simulation. The trajectory spirals away from the 
endemic state, with larger and larger oscillations until it hits the extinct state. The 
agreement between the stochastically simulated optimal path to extinction and the 
predicted optimal path is excellent. 

The curve of figure [T] is obtained by recasting the stochastic problem in a deter- 
ministic form. The evolution of the probability of finding a stochastic system in a 



given state X at a given time t is described by the master equation (van Kampen 



2007). Solving the master equation would provide a complete description of the 
time evolution of the stochastic system, but in general it is a difficult task to obtain 
explicit solutions for the master equation. Thus, one generally resorts to approxi- 
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mations to the solution; i.e., one considers an ansatz for the probability density. In 
this case, since extinction of finite populations is a rare event, we will be interested 
in examining events that occur in the tail of the probability distribution. Therefore, 
the distribution is assumed to take the form 

p(X,t) w exp(-NS(x)), (2.1) 

where p(X, t) is the probability density of finding the system in the state X at time 
t, N is the size of the population, x = X/iV is the normalized state (e.g., in an 
epidemic model, the fraction of the population in the various compartments), and 
S is a deterministic state function known in classical physics as the action. Equa- 



tion (2.1) describes the relationship between the action and the probability density 
and is based on an assumption for how the probability scales with the population 
size. The action is the negative of the natural log of the stationary probability dis- 
tribution divided by the population size. Therefore, the probability (if we normalize 
the population) is roughly given by the exponential of the action. Intuitively, equa- 



tion (2.1 1 expresses the assumption that the probability of occurrence of extreme 
events, such as extinction, lies in the tails of the probability distribution, which falls 
steeply away from the steady state. 

This approximation leads to a conserved quantity that is called the Hamilto- 



nian (Gang, 1987). From the Hamiltonian, one can find a set of conservative ordi- 
nary differential equations (ODEs) that are known as Hamilton's equations. These 
ODEs describe the time evolution of the system in terms of state variables x, which 
are associated with the population in each compartment. For the SIR example, x 
is the vector (S, I, R). In addition to the state variables, the equations contain con- 
jugate momenta variables, p x . The conjugate momenta, or noise, account for the 
uncertainty associated with the system being in a given state at a given time due to 
the stochastic interactions among the individuals of the population. These ODEs 
can be constructed from information in the master equation about the possible 
transitions and transition rates in the system. Details can be found in | Appendix A| 

Integration of the ODEs with the appropriate boundary conditions will then 
give the optimal evolution of the system under the influence of the noise. Boundary 
conditions are chosen to be fixed points of the system. A typical case is shown 
schematically in figure [2^,. Deterministically, the endemic state is attracting and 
the extinct state repelling. However, introducing stochasticity allows the system to 
leave the deterministic manifold along an unstable direction of the endemic state, 
corresponding to nonzero noise. Stochasticity leads to an additional extinct state 
which arises due to the general non-Gaussian nature of the noise. For the extinction 
process of figure[TJ boundary conditions were the system leaving the endemic steady 
state and asymptotically approaching the stochastic extinct state. 

In general, the optimal extinction path is an unstable dynamical object, and 
this reflects extinction as a rare event. This has led many authors to consider how 



extinction rates scale with respect to a parameter close to a bifurcation point (Do 



ering et al.\fM)S\ \K amenev fc Meerson[|2008l|Kamenev et aZ.l|2008l|Dykman et al 



2008), where the dynamics are very slow. For an epidemic model this means that 
the reproductive rate Rq should be greater than but very close to 1. However, most 
real diseases have Ro larger than 1.5, which translates into a faster growth rate 
from the extinct state. In general, in order to obtain analytic scaling results, one 
must obtain the ODEs for the optimal path either analytically (using the classical 
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theory of large fluctuations mentioned within this section and described in detail 
in Appendix A I or numerically (using shooting methods for boundary value prob- 



lems). This task may be impossible or extremely cumbersome, especially when the 
system is far from the bifurcation point. In the following section we demonstrate 
how to evolve naturally to the optimal path to extinction using a dynamical systems 
approach. 



3. Finite-time Lyapunov exponents 

We consider a velocity field which is defined over a finite time interval and is given by 
Hamilton's equations of motion. Such a velocity field, when starting from an initial 
condition, has a unique solution. The continuous dynamical system has quantities, 
known as Lyapunov exponents, which are associated with the trajectory of the 
system in an infinite time limit, and which measure the average growth rates of 
the linearized dynamics about the trajectory. To find the finite-time Lyapunov 
exponents (FTLE), one computes the Lyapunov exponents on a restricted finite 
time interval. For each initial condition, the exponents provide a measure of its 
sensitivity to small perturbations. Therefore, the FTLE is a measure of the local 
sensitivity to initial data. Details regarding the FTLE can be found in Sec. 2 of the 
electronic online supplement. 

The FTLE measurements can be shown to exhibit "ridges" of local maxima. The 
ridges of the field indicate the location of attracting (backward time FTLE field) 
and repelling (forward time FTLE field) structures. In two-dimensional (2D) space, 
the ridge is a curve which locally maximizes the FTLE field so that transverse to 
the ridge one finds the FTLE to be a local maximum. What is remarkable is that 
the FTLE ridges correspond to the optimal path trajectories, which we heuristically 
argue in Sec. 3 of the electronic online supplement. The basic idea is that since the 
optimal path is inherently unstable, the FTLE shows that, locally, the path is also 
the most sensitive to initial data. Figure [2]d shows a schematic that demonstrates 
why the optimal path has a local maximum to sensitivity. If one chooses an initial 
point on either side of the path near the endemic state, the two trajectories will 
separate exponentially in time. This is due to the fact that both extinct and endemic 
states are unstable, and the connecting trajectory defining the path is unstable as 
well. Any initial points starting near the optimal path will leave the neighborhood 
in short time. 



4. Evolving towards the optimal path using FTLE 

We now apply our theory of dynamical sensitivity to the problem of locating op- 
timal paths to extinction for several examples. We consider the case of internal 
fluctuations, where the noise is not known a priori, as well as the case of external 
noise. In each case, the interaction of the noise and state of the system begins by 
finding the equations of motion that describe the unstable flow. These equations 
of motion are then used to compute the ridges corresponding to maximum FTLE, 



which in turn correspond to the optimal extinction paths (Forgoston et al. 20101 
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(a) Extinction in a Branching - Annihilation Process 

For an example of a system with internal fluctuations which has an analytical 
solution, consider extinction in the stochastic branching-annihilation process 



A\2A and 2A 



(4.1) 



where A and fi > are constant reaction rates (Elgart & Kamenev 2004 As 



saf et al. 2008). Equation (4.1 1 is a single species birth-death process and can 



be thought of as a simplified form of the Verhulst logistic model for population 



growth (Nasell 2001 1. The mean field equation for the average number of individu- 
als n in the infinite population limit is given by h = Xn—jin 2 . The stochastic process 



given by equation (4.1) contains intrinsic noise which arises from the randomness 



of the reactions and the fact that the population consists of discrete individuals. 
This intrinsic noise can generate a rare sequence of events that causes the system to 
evolve to the extinct state. The probability P n (t) to observe, at time t, n individuals 
is governed by the master equation 

P n = | [(n + 2)(n + l)P n +2 - n(n - 1)P„] + X[(n- l)P n ^ - nP n ] . (4.2) 



The Hamiltonian associated with this system is 



H(q,p) = [X(l+p)-^(2+p)q qp 



(4.3) 



where q is a conjugate coordinate related to n through a transformation (Assaf 



given by 



et al. 2008|), and p plays the role of the momentum. The equations of motion are 

dH 



q=-g- =g[A(l + 2p)-M(l+p)g], 



P - 



dH 
dq 



(4.4a) 

p\ i i{2+p)q-X{l+p)]. (4.4b) 



The mean field is retrieved in equation (4.4a) when p = (no fluctuations 



or noise). The Hamiltonian has three zero-energy curves. The first is the mean- 
field zero-energy line p = (no fluctuations), which contains two unstable points 
hi = (p,q) = (0, A//i) and ho — [p,q) — (0,0). The second is the extinction line 
q = (trivial solution), which contains another unstable point h 2 — (p, q) = (—1,0). 
The third zero-energy curve is non-trivial and is 



q = 



2X(l+p) 
p{2+p) ■ 



(4.5) 



The segment of the curve given by equation ( |4.5[ ) which lies between — 1 < p < 
corresponds to a (heteroclinic) trajectory which exits, at t = — oo, the point hi 
along its unstable manifold and enters, at t = oo, the point h 2 along its stable 
manifold. This trajectory is the optimal path to extinction and describes the most 
probable sequence of events which evolves the system from a quasi-stationary state 



to extinction (Assaf et al. 2008). 
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To show that the FTLE evolves to the optimal path, we calculate the FTLE 



field using the system of Hamilton's equations given by equations (4.4a)-(4.4b). 
Figure [3^, shows both the forward and backward FTLE plot computed for the finite 
time T = 6, with A = 2.0 and /i = 0.5. In this example, as well as the following two 
examples, T was chosen to be sufficiently large so that one obtains a measurable 
exponential separation of trajectories. In figure [3ji, one can see that the optimal 
path to extinction is given by the ridge associated with the maximum FTLE. In 
fact, by overlaying the forward and backward FTLE fields, one can see all three 
zero-energy curves including the optimal path to extinction. Also shown in figure [3^, 
are the analytical solutions to the three zero-energy curves given by p = 0, q = 0, 



and equation (4.5). There is excellent agreement between the analytical solutions 



of all three curves and the ridges which are found through numerical computation 
of the FTLE flow fields. 

It is possible to compute analytically the action along the optimal path for a 



range of X//J, values. Using equation (4.5), it is easy to show that the action 5 is 



S = 2(1 -ln2)-. 



(4.6) 



It is clear from equation (4.6) that the action scales linearly with A//x. 



(b) SIS Epidemic Model - External Fluctuations 

We now consider the well-known problem of extinction in a Susceptible-Infectious- 
Susceptible (SIS) epidemiological model, which is a core model for the basis of many 
recurrent epidemics. The SIS model is described by the following system of equa- 
tions: 



S =fi — /j5 + 7/ — PIS, 

i = -(jt + i)i + pis, 



(4.7a) 
(4.7b) 



where pL denotes a constant birth and death rate, /3 represents the contact rate, and 
7 denotes the recovery rate. Assuming the total population size is constant and can 
be normalized to S + I = 1, then equations ( 4.7a )-( 4.7b ) can be rewritten as the 



the following one-dimensional (ID) process for the fraction of infectious individuals 
in the population: 

I = -(}i + i)I + pi(\-I). (4.8) 



The stochastic version of equation (4.8) is given as 



-(jt + 7)/ + /?/(! - 1) + a V (t) = F(I) + ar)(t), 



(4.9) 



where n(t) is uncorrelated Gaussian noise with zero mean and a is the standard 
deviation of the noise intensity. The noise models random migration to and from 



another population (Alonso et al.\ 2006 Doering et al. 2005) 



1 = (disease-free state) 
(/i + 7)//3 (endemic state). Using the Euler-Lagrange equation of 



Equation ( |4.9[ ) has two equilibrium points given by 



and 1=1 

motion (Goldstein et al. |2001 ) from the Lagrangian determined by equation (4.9) 
(L(I,I) = [ri(t)]' 2 = [I — F(I)] 2 ) along with the boundary conditions given by the 
extinct and endemic states, one finds that the optimal path to extinction (as well 
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as its counterpart path from the disease-free state to the endemic state) is given by 

i = ±F(i). 

As in the first example, one can numerically compute the optimal path to ex- 
tinction using the FTLE. Figure [3Jd shows the forward and backward FTLE plot 
computed for T = 5, with j3 = 5.0 and K = fi + 7 = 1.0. Note that we can consider 
the combination /1+7 since the Lagrangian depends only on the combination rather 
than on either parameter separately. In figure[3j}, one can see that the optimal path 
from the endemic state to the disease- free state is given by the ridge associated with 
the local maximum FTLE. Also shown in figure [3]d is the counterpart optimal path 
from the disease-free state to the endemic state (found by computing the backward 
FTLE field). In addition, the agreement with the analytical prediction is excellent, 
as shown in figure [3b. 



If one solves equation (4.9) for I(t) and substitutes both I(t) and I(t) into the 
Lagrangian, then one can analytically find an expression for the action along the 
optimal path. The expression for the action is a function of K and the reproductive 
number Rq and is given by 

where Rq = /3/k. 

(c) SIS Epidemic Model - Internal Fluctuations 

We next consider the ID stochastic version of the SIS epidemic model for a finite 
population with internal fluctuations using the transition rates found in Sec. 4 of 



the electronic online supplement. Using the formalism of Gang ( 19871, one then has 
the following Hamiltonian associated with this model: 

H(I,p) = ( M + l)I(e- p - 1) + 01(1 - I)(eP - 1), (4.11) 

where / is the fraction of infectious individuals and p is the momentum. Internal 
fluctuations arise from the random interactions of the population. Although there 
is no analytical solution for the optimal path to extinction, we can once again 
determine the optimal path by computing the FTLE flow field associated with 
this system. Figure [4] shows the forward FTLE plot computed using Hamilton's 
equation of motion for T = 10, with j3 — 2.0 and k = (1 + 7 = 1.0, and as in 
previous examples, the optimal path to extinction from the endemic state to the 
disease-free state is apparent. Note that the non-zero momentum corresponding 
to the extinct state qualitatively agrees with similar boundary conditions found 
in Dykman et al. ( 2008 ) and is associated with non-zero probability flux. 

Once again, it is possible to compute the action along the optimal path for a 
range of values of the reproductive number Rq. In contrast to the prior two exam- 
ples, here the action must be computed numerically. Moreover, even the optimal 
path must be found numerically using the FTLE plot generated for each value of 
R Q . 

Given an Rq, we computed the forward FTLE flow field using a grid resolution 
of 0.005 in both position and momentum space. To find the optimal path corre- 
sponding to the ridge of maximal FTLE values, we let the deterministic, endemic 
steady state be the starting point Zq of the path. Then a second point zi on the path 
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was found by locating the maximum of the FTLE values in an arc of radius 20 grid 
points spanning nearly tt radians for negative momentum values and originating at 
Zq. Subsequent points z„ + i were found by taking the maximum FTLE value on an 
arc of radius 10 grid points spanning nearly tt radians originating at the most recent 
point z n and centered around the vector z„ — z„_i until the extinction line was 
reached. The complete optimal path was estimated from the z„ using cubic spline 
interpolation. Once the optimal path had been found, the action was computed by 
numerically integrating the Lagrangian along this path. 

Figure [5^i shows the numerically computed action versus reproductive number 
over the range 1.1 < Rq < 20. The inset of figure [5^, shows the portion of figure [5^, 
from 1.1 < Rq < 2. Also shown in the inset of figure [5^, is an analytical, asymptotic 
scaling result that is valid for values of i?o close to unity. As can be see in figure 
there is good agreement between the numerically computed action and the ana- 
lytically computed action near Rq — 1. Details of the derivation of the analytical 
scaling law can be found in Sec. 5 of the electronic online supplement. 

Figure [5]d shows the numerically simulated mean extinction time versus re- 
productive number for the ID stochastic SIS model for a finite population (see 
) . Also shown in figure [HJd is the analytical prediction found using the previ- 
ously mentioned scaling law derived for values of Rq close to unity. As one can see, 
the agreement is excellent. 



5. Conclusions 



In all of the examples of Sec. |4j we have shown that the optimal path to extinc- 
tion is equated with having a (locally) maximal sensitivity to initial conditions. 
Even though there exist many possible paths to extinction, the dynamical systems 
approach converges to the path that maximizes extinction. The parameter values 
chosen for the three examples are such that the extinct and endemic states are far 
away from each other. Therefore, in general, there will be no possible approximate 



analytical treatment as performed in Dykman et at (2008). In addition, we have 



shown how to constructively compute numerically the action for a wide range of 
reproductive numbers. Our method allows for the computation of extinction times 
and can be extended to high-dimensional problems. 

Because the method is general, and unifies dynamical systems theory with the 
probability of extinction, we expect that any system found in other fields can be 
understood using this approach. Indeed, in problems of general extinction, it is 
now possible to evolve naturally to the optimal path, and thus predict the path 
that maximizes the probability to extinction. Future work in this area will include 
improved sampling methods to find the optimal path more efficiently in higher 
dimensional models. Specific applications of optimal path location in the future 
will include spatio-temporal extinction processes such as those that occur in prc- 



extinction in diseases such as dengue ( Cummings et al. 2005 ) 



vaccine measles (Keeling & Grenfell 1997 Finkenstadt et al. 2002) and multi-strain 



Finally, the optimal path method may lead to novel monitoring and control 
strategies. In one biological application, knowledge of the most probable path to 
extinction may help with monitoring an environment and with providing an es- 
timate of the likelihood of extinction based on where the population lies relative 
to the path. It is known that once a trajectory has left the neighborhood of the 
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endemic state, most paths to extinction occur near the optimal path. This phe- 
nomenon can be seen in figure [l] where the optimal path lies on the peak of the 
probability density of extinction prehistory. Therefore, the optimal path provides 
a good location to monitor the system for possible extinction events. Furthermore, 
although the momenta (noise effects) are not directly observable, it may be possible 
to infer them dynamically (Luchinsky et al. 2008) from time series data of observ- 



able quantities in conjunction with the equations for the time evolution of position 
and momentum variables. Knowledge of where a system lies in position-momentum 
space could provide an estimate for how quickly a desired epidemic extinction could 
occur or could provide the risk of extinction for a species one wishes to conserve. 

In yet another application, knowledge of the optimal path to extinction has the 
advantage that, in the presence of noise (that is estimated from data) and a known 
population of infectious individuals, it may be possible to develop better vaccine 
controls that reduce the time to extinction. Figure [2ja) shows a schematic of the 
path to extinction, where the extinct state is a saddle with stable and unstable 
directions. For many epidemic models, the extinct state can be shown to have a 
similar geometry of stable and unstable directions. An approach to the extinct state 
on the optimal path will lead to the fastest time to extinction. Moreover, since the 
extinct state has a saddle structure in the presence of noise, it may be controlled 



probabilistic techniques (Schwartz et al. 20041 



with projection methods (Schwartz & Triandaf 1994 Schwartz et al. 1997) or 



One can consider instituting a method of parameter control, where the param- 
eters could be vaccine levels or treatment of infectious individuals. Combined with 
the above-mentioned monitoring techniques, the control method will allow one to 
move an existing state that deviates from the optimal path closer to the optimal 
path as time evolves. By adjusting the parameters, we may target the stable direc- 



tions of the path when we are close to epidemic extinction (Schwartz & Triandaf 



1994 ) . Comparing observations with model predictions of the optimal path allow 



us to use controls to minimize the time to extinction. 
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Appendix A. Theory of large fluctuations 

Letting TV denote the population size, the state variables X £ M. n of the system 
describe the components of a population, while the random state transitions which 
govern the dynamics are described by the transition rates W(X,r), where r £ R n 
is an increment in the change of X. In the large population limit without any 
fluctuations, the mean field equations are given by the system X — ^ r rW(X, r). 

Consideration of the net change in increments of the state of the system results 
in the following master equation for the probability density p(X,t) of finding the 
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system in state X at time t: 

= J2 MX r; r)p(X -r,t)- W(X; r)p(X, t)] . (A 1) 

r 

The solution to the master equation (|A 1[) has a peak around a stable steady 



state in the limit of large population (N ^> 1) with width oc A 1 / 2 (Dykman et al. 



2008 Schwartz et al.. 20091. However, since we are interested in the probability of 



extinction, we will consider the tail of the distribution, which gives the probability 
of having a small number of individuals. The tail of the distribution can be obtained 



by employing the ansatz given by equation (2.1), which is an assumption for how 



the probability density scales with population size (Kubo et al. 1973 Dykman 
1990 Dykman et al. 2008). Equation (2.1) also implies that a maximum in the 



extinction probability can be found minimizing the action over a set of extinction 
paths starting from the stable endemic state. The assumption of this functional 
form for p allows the action to be derived from properties of the master equation. 

The density, p(X,t), can be found by substituting the ansatz given by equa- 
tion (2.1) into the master equation. The resulting equations for the action are 



given by the Hamilton- Jacobi equation for a Hamiltonian, H, given by dS/dt + 
H(x,dS/dx;t) — 0, where we have normalized the population components and 
rates respectively as x = X/N, w(x, r) = W(X, r)/N . 

Following the classical mechanics convention, define a conjugate momentum to 
the state space, x, by letting p — dS/dx and where H(x,p; t) is the classical Hamil- 



tonian (Kubo et al. 1973). The Hamiltonian, H, depends both on the state of the 



system, x, as well as the momentum, p, which provides an effective force due to 
stochastic fluctuations on the system. The Hamiltonian equations of motion provide 
trajectories in time of x(t) and p(t), and as such, describe a set of paths going from 
an initial point at time ti to some final point at time tf in (x,p) space. For a given 
path, the action is given by S = f. f p(t)x(t)alt, and as such determines the expo- 
nent of the probability of observing that path. (It should be noted that instead of 
using the Hamiltonian representation, one could use the Lagrangian representation 
L(x, x; t) — —H(x, p;t) + x ■ p, which results in an equivalent solution.) 

The action reveals much information about the probability evolution of the 
system, from scaling near bifurcation points in non-Gaussian processes to rates of 



extinction as a function of epidemiological parameters (Dykman 1990 Dykman 



et al. 2008). As already stated, in order to maximize the probability of extinction, 



one must minimize the action S. The minimizing formulation entails finding the 
solution to the Hamilton-Jacobi equation, which means that one must solve the 
2n-dimensional system of Hamilton's equations [x = d p H(x,p), p = —d x H(x,p)] 
for x and p, where the Hamiltonian is explicitly given as 



H(x,p, t) = w(x, r)[exp (pr) — 1] 



(A 2) 



The appropriate boundary conditions of the system are such that a solution starts 
at a non-zero state, such as an endemic state, and asymptotically approaches one 
or more zero components of the state vector, representing a disease free state. 
Therefore, a trajectory that is a solution to the two-point boundary value problem 
determines a path, which in turn yields the probability of going from the initial 
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state to the final state. The optimal path to extinction is the path which minimizes 
the action in either the Hamiltonian or Lagrangian representation. 

We compute the trajectory satisfying the Hamiltonian system that has as its 
asymptotic limits in time the endemic state as t — > — oo and the extinct state 
as t — > +00. The momentum p represents the force of the fluctuations on the 
population, and this momentum changes the stability of the equilibrium points. 
Both the endemic and extinct states have attracting and repelling directions for 
p ^= 0, as shown schematically in Figure [2] 

Given optimal path traj ectories (x(t),p(t)), th e action with correct limits is 



found by S(x) — J^pxdt (Goldstein et al. 2001) 
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Susceptibles 

Figure 1. Probability density of extinction prehistory and the optimal path to extinction 
for the stochastic SIR epidemiological model. Colors indicate the probability density (the 
colorbar values have been normalized, with lighter colors corresponding to higher proba- 
bility) for 20, 000 stochastic realizations. The results were computed using Monte Carlo 
simulations, and details of the sampling of the trajectories are described in the text. The 
curve is the numerically predicted optimal path to extinction. Note that the optimal path 
to extinction lies on the peak of the probability density of extinction prehistory. The pop- 
ulation is 3 x 10 6 individuals, with Ro « 15 (contact rate /3 = 1500, recovery rate 7 = 100, 
birth-death rate fj, — 0.2). 



population .v 




Figure 2. a) Schematic showing the fixed points and heteroclinic trajectories (trajectories 
connecting fixed points). Coordinate axes are dashed lines. The noise coordinate is indi- 
cated by the momentum (p) coordinate. The deterministic manifold (p = 0) is indicated in 
blue. Deterministically, the extinct state is repelling and endemic state is attracting. How- 
ever, the endemic state has unstable directions for nonzero noise (p 7^ 0), and the optimal 
path (red) is the trajectory carrying the system from the endemic state to a stochastic 
extinct state, b) Schematic showing the path from the endemic state (blue) to the extinct 
state (red). The optimal path leaves the endemic point along an unstable manifold, and 
connects to the extinct state along its stable manifold. The magenta and green dashed 
lines represent trajectories initially separated by the optimal path. The initial starting 
distance between trajectories near the endemic state expands exponentially in forward 
time (shown by the cyan lines). Locally, this shows that the finite-time Lyapunov measure 
of sensitivity with respect to initial data is maximal along the optimal path. 
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Figure 3. (a) Forward and backward FTLE flow fields for the branching-annihilation 

process computed using the Hamiltonian given by equation ( |4.3| l with A = 2.0 and pL = 0.5. 

The integration time is T = 6 with an integration step size of t = 0.1 and a grid resolution 

of 0.01 in both q and p (momentum). The three zero-energy curves are given by the ridges 

of maximal FTLE and are overlaid with the analytical solution of these curves given by 

p = 0, q — 0, and equation (4.5|. (b) Forward and backward FTLE flow fields for the 

SIS epidemic model with external fluctuations. The flow fields were computed using the 

Lagrangian given by equation ( |4.9[ ) with f3 — 5.0 and k = p + 7 = 1.0. The integration 

time is T — 5 with an integration step size of t — 0.1 and a grid resolution of 0.01 in 

both I and p (momentum). The optimal path to extinction from the endemic state to 

the disease-free state and its counterpart optimal path from the disease-free state to the 

endemic state (found by computing the backward FTLE field) are given by the ridges of 

maximal FTLE and are overlaid with the analytical solution of the optimal paths. 
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Figure 4. FTLE flow field for the SIS epidemic model with internal fluctuations computed 
using the Hamiltonian given by equation |4 . 1 1 1 with /3 = 2.0 and k — p, + 7 = 1.0. The 
integration time is T = 10 with an integration step size of t — 0.1 and a grid resolution of 
0.005 in both I and p (momentum). The optimal path to extinction is given by the ridge 
of maximal FTLE. 



Article submitted to Royal Society 



18 



LB. Schwartz, E. Forgoston, S. Bianco, and L.B. Shaw 




Figure 5. (a) Numerically computed action (integrated along the optimal path found 
numerically from the FTLE flow field) versus reproductive number Ro for the SIS epi- 
demic model with internal fluctuations. The inset shows a portion of the graph near 
Ro = 1. The numerically computed action is given by the black points, while the dashed, 
red curve shows an asymptotic scaling result that is valid near Ro = 1, and is given 
by S(Ro) = (Ro — l) 2 /[^o(l + Ro)] (see Sec. 5 of the electronic online supplement), 
(b) Numerically simulated (solid curve with black points) mean extinction time versus 
reproductive number for the SIS epidemic model with internal noise and the analytical 
prediction (dashed, red curve) found using the asymptotic scaling law that is valid near 
Ro = 1. 
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